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ABSTRACT 

Based on our newly developed methods and the XMM-Newton large program of SN1006, we extract 
and analyze the spectra from 3596 tessellated regions of this SNR each with 0.3-8 keV counts > 10^. 

For the first time, we map out multiple physical parameters, such as the temperature (/cT), electron 
density (rie), ionization parameter (net), ionization age {tion)^ metal abundances, as well as the radio- 
to-X-ray slope (a) and cutoff frequency {i'cutoff) of fho synchrotron emission. We construct probability 
distribution functions of kT and net, and model them with several Gaussians, in order to characterize 
the average thermal and ionization states of such an extended source. We construct equivalent width 
(EW) maps based on continuum interpolation with the spectral model of each regions. We then 
compare the EW maps of O VII, O VIII, O VII KJ — Ne, Mg, Si XIII, Si XIV, and S lines 
constructed with this method to those constructed with linear interpolation. We further extract 
spectra from larger regions to confirm the features revealed by parameter and EW maps, which are 
often not directly detectable on X-ray intensity images. Eor example, O abundance is consistent with 
solar across the SNR, except for a low-abundance hole in the center. This “O Hole” has enhanced 
O VII KS — ( and Ee emissions, indicating recently reverse shocked ejecta, but also has the highest riet, 
indicating forward shocked ISM. Therefore, a multi-temperature model is needed to decompose these 
components. The asymmetric metal distributions suggest there is either an asymmetric explosion of 
the SN or an asymmetric distribution of the ISM. 

Subject headings: ISM: supernova remnants; acceleration of particles; shock waves; X-rays: ISM; 
methods: data analysis; (ISM:) cosmic rays. 

1. INTRODUCTION 

In young supernova remnants (SNRs), both shocked thermal plasm a emission and non-therma l emission fr om ac¬ 
celerated particles typically peak in the 0.5-10 keV X-ray band (e.g., iRnynold ^ 120081: IVinkI [20T1: 1 S land [201^ . The 
X-ray properties of these emission components often show clear spatial variations. Therefore, spatially resolved X-ray 
observations of young SNRs play important roles in our understanding of many scientific issues, such as the explosion 
mechanisms and progenitors of supernova (SN), the microphysics involved in particle acceleration and magnetic field 
amplification, the heating of electrons an d ions at the shock, as well as the abundances and the distribution of fresh 
nucleosynthesis products (see IVinljl2Q12l for a recent review). 

Current facilities (in particular. X-ray CCDs) on board space X-ray telescopes such as Chandra and XMM-Newton 
have the ability to simultaneously record the spatial and spectral information of the collected photons. They could thus 
be used for spatially resolved spectroscopy analysis to study the spatial distribution of X-ray properties across young 
SNRs. Howe ver, methods com monly used in X-ray data analysis, such as spectral analysis of i ndividual interestin g 
regions (e.g., iChen et al.ll2QQ8l ) and equivalent width (EW) map of strong emission lines (e.g., iHwang et al.ll2QQQl ). 
usually do not make full use of the information contained in the X-ray CCD data. What we want, and are often 
contained in the X-ray data, are the spatial distributions of various physical properties of the thermal and non-thermal 
X-ray emission across the SNRs. 

In this paper, we will introduce new techniques to conduct spatially resolved spectroscopy analysis. The basic idea is 
to map out the spectral parameters in small tessellated regions which individually contain enough photons for spectral 
analy sis. Such techniques have initially been developed for optical integral field observations (l^ppellari CopinI 
l2QQ3f ) , and have later been applied to X-ray observations (jDiehl X Statlerl I2QQ6I: iBroos et~an I2Q1QD , especially ^ 
resoly ing interesting features in t he temperature/metallicity maps of massiye galaxy clusters (e.g., iRandall et al.l 
l2QQ8f) or interacting galaxies (e.g., iHodges-Klnck fc RnynoldsIl2Q12[ ) . SNRs t ypically hav e more comp licated X-ray 
spectra which often include ionization non-equilibr ium thermal plasm a (e.g., IVinkI |2QT^ iSland 1^141 ) and/or non- 
thermal emission with varying spectral shapes (e.g.. iMiceli et alll2Q14f ). It is therefore more difficult to automatically 
decompose different spectral components and model the spectra extracted from a large number of tessellated regions. 
Eurthermore, existing sp atially resolved spectroscopy analysis of SNRs of ten adopt equal-sized (e.g..lLn fc Aschenba^ 
l2QQQl : iLopez et al.l l2Q13h or adaptively-binned box-shaped meshes (e.g., iCassam-Chenai et al.ll2QQ3 ). which are less 
efficient in resolving fine structures on physical parameter maps. 

The remnant of the supernova AD1006 (SN1006) is one of the few historical SNRs with human records of its exact 
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birth date (e.g., iStephens^ 1201 Q[) . It is widely accepted that SN1006 is the remnant of a Type la supernova event, 
based on its high Galactic latitude location \ b = 14.6°, apparently isolated from any star formation re gions], the lack 


of any visible central compact sources (e.g., iPve et all 119811: iJones fc Pvelll989l : iBnrleig! 


of iro n absorption lines from the ultraviolet (UV) spectra of background sources (e.g. 


1 et al.ll2QQQ), the detection 


119931 ). the metal abundances inferred from the soft X-ray emission lines (iKovama et a 


Fesen et al.l 119881: IWu et al.l 

] [2(M lU^ida et al.l IMl) . 


2Q1Q[ ). and the historical 


the lack of time variation due to the clumpy of th e surrounding me dium (jKatsuda et al. 
records that it re mains visible for several years (e .g. jStephensonll2QlQD. Proper m otion measurements in optical (e.g. 


/inkier et al.N2Q14f ). combined with the shock velocity 


iLong et al]ll988f ). radio (e.g.. iMoffett et aTlll993D . or X-ray ^ _ 

as in f erred from the measurements of the width and ratio of some optical/UV emission lines (e.g., iKirshner et al.l 
119871 : iLaming et al.l Il9^ iGhavamian et al.ll2QQ2D. or the expan ding velocity inferred from the broadening of some 
UV absorption lines of background sources (e.g., IWu et al1ll993[ ). place this remnant at a distance of 2.18 T 0. 08 kpc 
(| Winkler et al.l [200^ . well consistent with the distance obtained from H I observations (jPubner et al.l[2QQl ). The 
apparent radio/X-ray diameter of SN1006 is ^ 30', or ^ 19 pc at this distance. 

Similar as other young SNRs, the multi-wavelength properties of SN1006 also show clear spatial variations. First, 
SN1006 ca n be apparently divided into two distinguisha ble parts. The nor theast (NE) and southwest (SW) lobes ar e 
radio fe.g .. iRevnolds fc Gilmorel[l986l: iDver et al.ll2QQ9[ ). hard X-ray (e.g., iKovama et al.l [19^ IWinkler et al.ll2Qf^ . 
and TeV (lAcero et al.ll2QlQD br ight, and the axis connecting the bright non-thermal regions is roughly aligned with the 
Galactic plane (jGaenslerlll998[ ). They are thought to be do minated by synchrotron emission of relativistic electrons 
accelerated at the SNR blast wave (e.g.. IKovama et 'allll995f ). On the other hand, most of the SNR interior, including 
the n orthwest (NW) shell , are dominated by thermal emission characterized by strong emission lines of heavy elements 
(e.g.. IKovama et aIll2QQ8l:IUchida et al.ll2Q13f ). Second, not only the relative contribution of the two major components, 
but also the spectral properties of the non-thermal emission within the NE/SW lobes and the thermal emission within 
the SNR interio r, show significant variations. These variations indicate the spatial vyiation of particle acceleration 
conditions (e.g., lR.othenflug et al.l[2004 iCassam-Chenai et al.ll2QQ8|: [X iceli et al.lI2Q13LI2Q14D. as well as the thermal 
and chemical properties across the entire remnant (e.g.. lUchida et alT 20131: IWinkler et al.l2Ql^ . which await better 
characterizations. Third, the multi-wavelength properties of SN1006 also sh ow noticeable spatial variations. In addition 
to the two promin ent non-therma l lobes, the spatial v ariation of Ha (e.g., IWinkler et al. I2QQ3I: lRaymonc~eF~aIll2QQ7|: 
Nikolic et al.ll2Ql^ . infra-red (IR: IWinkler et al]l2Q13[ ). and H I 21-cm line emissions (jPubner et al.ll2QQ2l: iMiceli et al l 


201 4 ) on both large (comparable to the size of the SNR) and small (comparable to the size of some prominent 


features) scales further indicates the considerable spatial variation of the density of the ambient interstellar medium 
(ISM), although SN1006 is believed to evolve in a relatively low-density and uniform environment. 

We herein study SN1006 with our newly developed spatially resolved spectroscopy analysis method. We have 
obtained high quality X-ray data of SN1006 through our XMM-Newton Large Program (LP) and some archival 
observations. The high sensitivity and large field of view (FOV) of XMM-Newton help us to collect enough photons 
over the entire remnant, thus make SN1006 the best candidate to test how the spatially resolved spectroscopy analysis 
with our new techniques could improve our understanding of young SNRs. The paper is organized as follows: In ^ 
we briefly introduce our XMM-Newton LP and also the archival data used in this project. Basic data calibration is 
also detailed in this section. In ^ we describe our new methods developed to conduct spatially resolved spectroscopy 
analysis. The parameter maps and other products of this analysis are presented and further discussed in 0 We 
summarize the main results and conclusions in ^ More discussions on the thermal and non-thermal emissions of this 
remnant will be presented in companion papers. 


2. OBSERVATIONS AND DATA CALIBRATION 
2.1. XMM-Newton Large Program and Archival Data 

In this work, we analyze the data obtained from the XMM-Newton LP of SN1006 (PI: A. Decourchelle, ^ 700 ks 
of total exposure time). These observations were all performed with the “Medium” filter by using the full-frame 
(FF) mode for the MOS (Metal Oxide Semi-conductor) cameras and the extended full-frame (EFF) mode for the 
PN cameras. We also select archival XMM-Newton EPIG (the European Photon Imaging Gamera) observations with 
pointing positions within 30' from the center of SN1006. Only observations with at least one of the EPIG cameras 
operating in either FF or EFF mode are considered. All these observations are also carried out with the “Medium” 
filter. Information of all the selected o bservations is summarized in Tab l e [H S ome analysis of these LP data have 
already been published in IMiceli et al.l (j2012L I2013L l2014l ): iBroersen et al.l (j2013f ). 


2.2. Data calibration 

We reduce the data based on XMM-Newton Science Analysis Software (SAS) vl2.0.1. For each observation listed in 
Table [H the Observation Data Files (ODF) for each of the EPIC instruments (MOS-1, MOS-2 and PN) are reprocessed 
using the SAS tasks emchain and epproc. We identify and tag the low-energy noise in the MOS CCDs in anomalous 
states, i.e., with an elevated event rat es between 0 and 1 keV, using the SAS task emtaglenoise, following the algorithm 
described in iKuntz fc SnowdenI (j2QQ8[ ). All events in all CCDs tagged as noisy are then filtered out. CCDs in anomalous 
states with this low-energy noise are also summarized in Table [TJ Each dataset is further screened for periods of soft 
proton flaring through the creation of full-held (but remove brightest point-like sources which may be time variable) 
light curves in broad band (0.3-12 keV for MOS and 0.3-14 keV for PN). Good Time Interval (GTI) selections are made 
by setting a threshold with 3 a clipping for each instrument. The resulting effective exposure times for each instrument 
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TABLE 1 

XMM-Newton observations of SN1006 used in this work 


ObsID 

Start Date 

noise CCD 


tliv,M2 

tliv,PN 

A//,Ml 

teff,M2 

teff,PN 

0077340101 

2001-08-10 

_ 

64862 

64910 

55980 

30716 

31148 

20896 

0077340201 

2001-08-10 

- 

57180 

57213 

44280 

24678 

24491 

18969 

0111090101 

2000-08-20 

- 

7425 

7430 

2954 

7325 

7430 

2954 

0111090301 

2000-08-17 

- 

5009 

5062 

899 

1678 

1679 

0 

0111090601 

2001-08-08 

- 

15770 

15776 

10387 

7560 

7957 

4337 

0143980201 

2003-08-14 

- 

30077 

30089 

23184 

16896 

16904 

11915 

0202590101 

2004-02-10 

- 

42796 

42886 

35155 

29421 

31563 

19733 

0306660101 

2005-08-21 

4 of MOS-1, 5 of MOS-2 

33424 

33446 

25526 

9546 

10528 

6582 

0555630101 

2008-08-22 

5 of MOS-2 

44851 

44884 

34437 

44551 

44781 

26517 

0555630201 

2008-08-04 

4 of MOS-1 

107848 

107885 

84401 

94477 

95999 

52548 

0555630301 

2009-02-19 

- 

119926 

119861 

93164 

89856 

91656 

59983 

0555630401 

2009-01-28 

4 of MOS-1 

102974 

103016 

81374 

75749 

81992 

41791 

0555630501 

2008-07-31 

5 of MOS-2 

125449 

125430 

100663 

85498 

96372 

51187 

0555631001 

2008-08-28 

5 of MOS-2 

64641 

64665 

50191 

60775 

61571 

44400 

0653860101 

2010-08-28 

4 and 5 of MOS-1 

124863 

124952 

101182 

104396 

105661 

77326 

Total 

- 

- 

947095 

947505 

743777 

683122 

709731 

439140 


Note. — The third column summarizes the identification number of MOS CCD(s) with low-energy noise. tuy^Mi, 
tiiv,M2, and tiiy^PN are the dead-time corrected on time (LIVETIME) for EPIC MOS-1, MOS-2, and PN, respectively. 
teff,Mi, teff,M2, and teff,PN are the effective exposure times of EPIC MOS-1, MOS-2, and PN after background ffare 
removal. The last row is the total exposure time of individual cameras of all the data used in this work. All the exposure 
times are in unit of second. 

after this background flare filtering are listed in Table [H The total effective exposure times of all the observations 
for MOS-1, MOS-2, and PN are 683, 710 and 439 ks, respectively. Out-of-Time (OoT) events have a non-negligible 
impact on the PN data and are also reprocessed with the SAS task epchain and filtered in an identical way to the 
primary PN datasets. This OoT events are further subtracted in the following imaging and spectral analyses. 

The instrument (telescope+detector) response is not fiat, i.e., the effective area at a given energy depends on the 
position in the focal plane. This vignetting effect is extremely important for the analysis of inhomogeneous extended 
source such as SN1006. We apply the SAS task evigweight to both the OoT and the primary events lists, weighting 
each EPIC events with inverse effective area over one exposure, so that the derived event list is equivalent to what 
one would get for a flat instrument. This allows the use of single on-axis Ancillary Response Files (ARE) for each 
instrument, and does not require further effective area corrections in the following imaging and spectral analyses. 

The source event maps and the instrument exposure maps from each observation and each camera are constructed 
in several bands. In order to correct for the OoT events in PN observations, we also construct OoT images and 
scale them to the expected OoT event fraction (6.3% in FF and 2.2% in EFF modes respectively, according to 
XMM-Newton Users Handbook: http://xmm.esac.esa. int/external/xmm_user_support/documeiitatioii/uhb/ 
index.html). These OoT images are then subtracted from the raw PN images. A subtraction of the instrumen¬ 
tal background is applied to images from all instruments by making use of the EPIC Filter Wheel Closed (FWC) 
data obtained from the XMM-Newton background analysis website (http: //xmm2. esac. esa. iiit/external/xmm_sw_ 
cal/background/f ilter_closed/index.shtml). These FWC data are reprocessed and renormalized to match the 
I0-I2/I2-I4 keV (for MOS/PN) counts number of the source images before subtraction. 

The source, exposure, and background maps from different observations are then combined together using the SAS 
task emosaic. All the images are binned to a pixel size of 3.2", slightly smaller than the FWHM (Full Width at 
Half Maximum) of the XMM-Newton mirror on-axis Point Spread Function (PSF; ^ 6"). Therefore, we do not lose 
resolution in constructing images. The background-subtracted event map is further adaptively smoothed with the SAS 
task asmooth to a desired signal-to-noise ratio of 5. The exposure map is smoothed according to the same template as 
the event map. We then produce final background-subtracted and exposure-corrected flux images with these smoothed 
images (e.g.. Fig. [T]). 

Because there exist a lot of soft X-ray knots in SNI006, which may be mis-identified as point sources in usual source 
detection tools, we conduct point-source detection only in the hard X-ray band (2-8 keV). We create a constant PSF 
map of 10" for the mosaiced images, which is good enough for a simple detection and removal of them for the study 
of diffuse emission. We then adopt this PSF map to a standard wave detection tool wavdeteet. We finally inspect the 
detected point-like sources visually to remove obvious false detections. The point sources are only removed in extracting 
the background spectra (Appendix [A|) . There are many faint point-like sources projected inside the SNR. The X-ray 
emission of most of these sources peak at hard X-ray so do not contribute significantly at < 2 keV which we most 
concern in this paper. Therefore, we do not remove these point sources from the source spectra. We caution that some 
peculiar features on the parameter maps shown in the following sections may be caused by foreground or background 
point sources instead of the emission from the SNR itself. This contamination is in general not significant, but could be 
important when the truly diffuse emission is faint (e.g., close to the two point sources at RA, Dec ~ I5^04’^20^, —42°02' 
outside the southeast rim of the SNR; Fig. [2j). 

3. SPATIALLY-RESOLVED SPECTROSCOPY 
3.1. Mapping the speetral parameters: general proeedure 
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Fig. 1.— Broad band tricolor images (Red: 0.3-1 keV; Green: 1-2 keV; Blue: 2-8 keV) of the entire field of view of the XMM-Newton 
observations of SN1006. The white solid contour shows the region to extract spectra for the whole remnant and to do spatially resolved 
spectroscopy analysis. The white dashed annulus show the region to extract the sky background spectra, and is divided into four quadrants 
to study the possible azimuthal variation of the sky background (Appendix |AJ . To clarify, removed point-like sources are not shown. 


X-ray observations typically provide us with multi-dimensional data, i.e., by recording the (x,y) position on the focal 
plane, the arrival time, and the energy of the photons. This multi-dimensional data allows us to study the spatial 
distribution of many physical parameters, in addition to simply constructing broad- or narrow-band intensity images. 
Spatially-resolved spectroscopy represent such techniques to map out the physical parameters of X-ray bright extended 
sources. Different from previous studies which often extract spectra from several interesting regions of an extended 
source, we herein introduce a new method to directly map out the spectral analysis parameters in r nany tessellated 
regions. Similar techniques have been discussed in several literatures. For example. iRandall et al.l (j2QQ8[) used two 
different methods (the over-sampling of the image at each pixel or the tessellated mesh) to map out the temperature 
distribution around the Virgo cluster galaxy M86. 

3.1.1. Creating tessellated meshes 

We adopt a new algorithm in order to construct tessellated meshes adaptively. We first find the brightest unbinned 
pixel in the original point-source-removed broad-band (0.3-8 keV) counts image and put this pixel into a new mesh. 
If the pixels already included in the mesh contain total (MOS-l+MOS-2-fPN) counts number above our threshold of 
10^ counts, we create one polygon region containing all the pixels in the mesh. If the pixels contain less counts, we 
then add the brightest pixel surrounding the previous added pixel to the mesh and compare the total counts to the 
threshold again. This process is repeated until the mesh contains more counts than the threshold or no surrounding 
pixel is unbinned. Isolated pixels between different meshes are leaving unbinned during the creation of the meshes. 
These unbinned pixels are finally individually merged into the neighboring meshes with the smallest average distance 
to them. This algorithm creates tessellated polygon regions each has at least a counts number of the input threshold, 
while each pixel only belongs to one single region (the meshes). An example of the created meshes (3596 tessellated 
regions), as will be adopted in the following analysis, is shown in Fig. [2l The smallest meshes have typical diameter 
of 4 pixel, or ~ 13", larger than the FWHM of the PSF. Therefore, the tail of the PSF should have little effect on the 
spectral analysis. 


3.1.2. Extraeting speetra of individual regions 

We extract spectra of MOS-1, MOS-2, and PN from each of the polygon regions of each observations containing 
positive counts numbers. We then weight these spectra as well as the corresponding background and response files 
and stack them together for MOS-1, MOS-2, and PN, respectively. Typically, the spectrum of each instrument from 
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Fig. 2.— Tessellated meshes used for spatially resolved spectroscopy analysis (the high resolution run) overlaid on the 0.3-8 keV image 
of SN1006. There are in total 3596 regions, each containing > 10"^ counts from the combination of MOS-1, MOS-2, and PN. The four black 
meshes are the regions used to extract the sample spectra in Fig. [S] with the region numbers denoted beside. 


a single region contains several thousand net counts, enough for the purpose of roughly characterizing the soft X-ray 
spectral energy distribution. To save computer time and space, we do not create response files [the Redistribution 
Matrix File (RMF) and the ARF] for all the spectra extracted for the three instruments of the 15 observations from the 
3596 regions 10^ spectra and the corresponding RMF; hereafter referred as the high resolution run). Instead, we 
create lower resolution meshes (with 182 regions in total, hereafter referred as the low resolution run), thus with higher 
counting statistic, and create RMFs for the spectra extracted from each regions of the three instruments of all the 
observations. The RMFs are then weighted according to the effective exposure time of individual spectra and merged 
together to create single weighted RMF for each region. These RMFs are taken as templates in the high resolution 
run. In the high resolution run, for each region, we find the nearest region from it in the low resolution run and directly 
use the existing RMF from this region. Because we have already corrected the vignetting in the calibration ( ^2.2p . we 
use identical on-axis ARF of each instrument for all the regions. In both low and high resolution runs, we still extract 
the source and background spectra separately for individual regions, instruments, and observations, and stack them 
according to their effective exposure time and area scale, similar as the standard way typically taken by many other 
authors. The use of template RMF may cause some biases in the spectral analysis, but we have double checked all 
the prominent features in the spectra extracted from the tess ellate d meshes, by examining the spectra extracted from 
larger regions with the response files individually generated ( ^4.ip . 

3.2. Spectral modeling of individual regions 
3.2.1. Spectral model 

Analysis of the sky background is presented in Appendix [A] We add the sky background components (with all the 
parameters fixed) to the model of the source spectra. In order to do this, the normalization of the three sky background 
components (VMEKAL, power law, and MEKAL) are rescaled according to the ratio of the effective areas of the source 
and sky background regions. In the analysis of both source and sky background spectra, we adopt the EWC data as 
the instrumental background and have subtracted them before spectral analysis. 

We fit the source spectra with a “VNEI+SRCUT” model. The VNEI component (Non-Equilibrium Ionization 
collisional plasma model, with parameters describing the relative abundances between different elements) describes 
the thermal plasma contribution, assuming a constant temperature and single ionization parameter. The SRCUT 
component describes the syn chrotron emission from an exponentially cutoff power-law distribution of electrons in a 
homogeneous magnetic field (|Rnvnolds fc KeohandflOQ^ . 

Both the VNEI and SRCUT components are subjected to foreground absorption described with a WABS model 
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(photo-electric absorption using Wisconsin cross-sections). It is possible that the distribution of the ambient ISM 
is inhomogenons , which could affect the absorption and thus the shape of the soft X-ray spectra. However, based 
on iDnbner et al.l (|2QQ2[ )’s H I column density map in the surrounding r egion of SN1006, the variation of absorbing 
H I column density (A^h) is at most ^ 25% (also see iMiceli et al.l 120141 ). Because A^h in the direction of SN1006 is 
quite low, such small variations on A^h does not significantly affect the spec tral fitting results. We therefore fix the 
absorption column density at the foreground value of A^h = 6.8 x 10 ^^ cm“^ (IDnbner et al.ll 2 QQ 2 li. _ 

For the VNEI component, we adopt the Atomic Data for Astrophysicists (AtomDB: iFoster et al.ll 2 Q 12 f ). We caution 
that the VNEI model only has a single ionization parameter (riet), so does not take account the possibly different 
ionization age of different elements or a distribution of ionization age in the post-shock region. The ionization age 
distribution of the post-shock plasma, and/or the different temperatures of the ions and electrons, are better described 
with a plane-parallel shock plasma model (YPSH OCK) or an X-ray spectral model describing SNRs in the Sedov- 
Taylor phase (Sedov model; iBorkowski et al l[ 2 QQlf ). However, spectral fitting with such models take too long computer 
time for the 3596 regions, so are not adopted here. 

Our single temperature (1-T) VNEI model may also fail to describe the possible multi-temperature structure of the 
plasma (e.g., lUchida et al.l[201^ . However, most of the spectra from individual meshes have poor counting statistic 
at high energy (e.g., at the Si line band), so have poor constraint on the properties of some possible high-temperature 
components. Therefore, the spectral model adopted in the present paper is only aimed at roughly decomposing the 
thermal and non-thermal contributions, as well as characterizing the average thermal and ionization states of the 
plasma. We will briefly discuss the thermal structure of the plasma in ^4.2.11 Eurther discussions on modeling the 
thermal spectra of SN1006 based on higher-quality spectra but lower-resolution meshes will be presented in companion 
papers. 

We set the O, Ne, Mg, Si, and Ee abundances free, the S abundance equal to Si, and Ni abundance equal to Ee. 
All these elements produce strong emission lines in soft X-ray, but Ee and Ni do not have well separated emission line 
bumps at < 2 keV, in which our spectra typically have enough counts (e.g., Eig. [3]). 

We further add into our models three narrow gaussian lines (line width fixed at 10“^ keV) at 0.714 keV, 0.723 keV, 
and 0.730 keV to account for the high er level transition s (K^, Ke, K^) of O VH lines which are not included in the 
VNEI code (e.g., Eig. [3K,d). Eollowing lYamagnchi et all (j2008f ). we assume Ke/K(5=KC/Ke=0.5, so there is only one 
free parameter of this component (the normalization of the O VH K(5 line). At temperature kT > 0.6 keV, the O VH 
lines are typically weak compared to the O VHI lines. At lower temperatures where the O VH line is dominant, the 
line flux decreases rapidly along the K-shell transition series, so in most cases only Ka — 7 transitions are important. 
However, SN1006 has a low ionization age and fairly high temperature It is thus possible that these high level 

transitions are fairly strong in some regions, the spectra of which show clear residual at 0.7-0.75 keV without adding 
these gaussian lines. Alternatively, the spectral fitting residual at 0.7-0.75 keV could also be explained as the Ee L 
shell transitions from low ionization ions (less than Ee^^+; e.g., IWarren fc Hugh^ 120041: lUchida et al.ll 20 f^ . We will 
further discuss this possibility in ^4.3.51 Nevertheless, for the convenience of presentation, we will call these added 
gaussian lines the O VH K(5 — lines throughout the paper, regardless of its real origin from oxygen or iron. 

In the SRCUT model, we fit the radio-to-X- ray photon index a and the cutoff frequency I'cutoffi but fix the 
normalization us ing the radio flux obtained fromlDyer et ahl (l2009l). We a ssume a constant radio spectral slope of 0.5 
{a ^ 0.45 — 0 . 6: iDver et al.l[ 2001 l: lAcero etldlDOOTl iKatsnda et al.ll 2010 [ ). in order to convert to the flux at 1.4 GHz 
in iDver et ^ (j2009l ) to the flux at 1 GHz adopted in the SRGUT model. Small variation of the assumed a does not 
significantly affect our spectral analysis results. 

The spectra of the three EPIG camera (MOS-1, MOS- 2 , and PN) are jointly fitted with all the parameters linked, 
except for a constant normalization factor subjected to all the model components. This normalization factor is used to 
account for the possible difference in area scale and calibration bias of the spectra extracted from different instruments. 
As shown in Eig. |4t^,b, it is close to 1 over the entire remnant, with slightly larger bias appears near the gaps between 
different GGD chips. This bias will not significantly affect our results, however. Therefore, in the following analyses 
and discussions of this paper, we will directly use the parameters from jointly fitting the MOS-1, MOS- 2 , and PN 
spectra, without considering the normalization factor. 

We also add a gain correction to the response file of the PN spectrum (adopting the “gai n” model in XSpec) . This 
is aiming at accounting for the deficiency in the low energy calibration of the PN camera (jPennerl et al.l [20041 ). The 
slope of “gain” is fixed at 1 and the offset is set free. 

Such a “VNEI+SRGUT-t-0 VH K(5, e,C line+Background” model typically gives acceptable spectral fitting results. 
Examples of spectra extracted from individual regions dominated by different spectral features are shown in Eig. [3] 
and the best-fit spectral parameters of them are summarized in Table [2] together with their statistical errors. The 
maximum y^/d.o.f for individual tessellated meshes is ^ 2.81, and > 93% of the regions have y^/d.o.f < 2.0 (Eig.Hh). 
We have individually inspected all the spectra with y^/d.o.f > 2.0, and refit the spectra if necessary. In most of the 
cases, the large y^/d.o.f value could be explained by the too simple 1-T spectral model adopted in the present paper. 

3.2.2. Derived parameters 

We further derive some physical parameters based on the parameters directly obtained from spectral fitting. We 
first estimate the electron number density (rie) from the emission measure of the plasma. In the VNEI model, the 

normalization is defined as: norm = f where Da is the distance to the source in cm, uh is the 

hydrogen number density, and dV is the volume of the projected spectral analysis regions. Adopting an outer shell 
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Energy (keV) Energy (keV) 


Fig. 3. — Examples of spectral fitting of individual tessellated regions (region numbers denoted in Fig. [2])- The black, red, and green 
data points are the instrumental background-subtracted spectra of MOS-1, MOS-2, and PN, while the solid curves with the corresponding 
colors are their best-fit models. The vertical axis of the upper half of each panels is the intensity per unit effective area (different from 
other spectra shown in this paper), so the spectra of different instruments have roughly identical intensities. The blue solid curves are the 
source components (VNEI, SRCUT, and Gaussian lines) of the PN spectrum, while the dashed curves are the sky background components 
(the local hot bubble. Galactic halo, and distant AGNs; Appendix Ia)) . The lower panels show the residuals in terms of a. Panel (a) is 
dominated by non-thermal emission, while panels (b-d) are dominated by thermal emission, with strong Ne, Si, and O lines, respectively. 
Panel (d) also has a relatively large contribution from the O VII K(5, e, ^ lines (the three Gaussian lines marked in blue solid curves), while 
this component is too weak to appear on the plot in panels (b) and (c). 


radius of 9.5 pc and an ambient medium density of no ^ 0.05 cm“^ (e.g.. lMiceli et al.ll2012[) and solar abundance, 
the total swept up ISM mass by SN1006 will be ^ 5 Mq, not much higher than the total ejecta mass of a Type la SN. 
Therefore, we may expect a significant metal enrichment by SN ejecta, as indicated by the significant spatial variations 
of different elements (see gj). However, in the present paper, we do not decompose the ISM and ejecta components in 
our simple 1-T thermal plasma model. In most of the spectra dominated by thermal emission, O lines are the most 
prominent (e.g.. Fig. 0), and the fitted O abundances are close to solar ( ^4.3.ip . Therefore, for the mixed ISM plus 

















































































































































































Fig. 4.— Panels (a) and (b): Normalization factor of MOS-2 (a) and PN (b) to the MOS-1 spectra, in order to account for the possible 
difference between the spectra extracted from different instruments. Panel (c): (d.o.f is the degree of freedom) of the spectral 

fitting of each tessellated regions. 


TABLE 2 

Examples of spectral fitting parameters and errors 


Region 

kT 

O 

Ne 

Mg 

Si 

Fe 

log net 

normvNEl 

F 

^cutoff 


keV 

solar 

solar 

solar 

solar 

solar 

log(cm“^s) 



10® GHz 

100 

> 0.22 

0.84tl:ll 

< 0.26 

< 34.4 

< 28.0 

< 13.5 

Q 40 + 1-22 
^•^^-0.33 

i.UU_o 49 

0.49 ±0.01 

46.411,®3 

1480 

n oq + 0.09 


< 7.43 

< 3.16 

< 0.17 

> 10.07 

34.21®;® 

< 0.15 

0.9ll°;™ 

3038 

1 47+0.33 
-0.22 

u.oy_o 26 

0.28 ±0.06 0.59to,il 

25.llfo" 

1 p:7+2.15 
-0.51 

Q ^7+0.05 
-0.08 

4 41+1.69 
^•^-'■-2.07 

< 0.14 

n 73+0.16 
'"^-O.OS 

3534 

4.ou_o 23 

< 19.5 

6.361N? 

50.7tlti 

420l®«o 

< 6.96 

Q 43+0.06 
^•^'^-0.04 

0'-^O-0.02 

0 1 in+0.032 
U.iiD_o 004 

0 «3+6-21 


Note. — Some key spectral fitting parameters and their 90% confidence errors of the spectra shown in Fig. normvNEi is in unit of 


4 [-D^ ( 1 + )]^ rieniidV. See the text in ^3.2.2| for details. 


ejecta spectra, we assume solar abundance in the calculation of Ue for simplicity. The assumption of solar abundance 
links the electron and hydrogen densities by Ue = 1 . 2 Inn- 

We adopt the spherical shell geometric model described in iMiceli et al.l (|2Q12[) to estimate the volume of the spectral 
regions. We assume SN1006 has an outer radius of 11 pc, slightly larger than the average X-ray radius of ^ 9.5 pc 
(2 ]) 5 ill order to cover all the X-ray emitting regions (the two non-thermal “ears” are slightly more extended). We 
also equalize each tessellated mesh to a sector-shaped region with the same area for the convenience of calculating its 
volume. There are two poorly constrained parameters of this model: the thickness of the thermal X-ray emitting shell 
and the volume filling factor / within the X-ray emitting regions. We assume different shell thickness and estimate the 
plasma density accordingly. In Fig. [5l we plot rig of some outermost regions dominated by thermal emission against 
the assumed shell thickness. 

The choice of regions dominated by thermal emission reduces the uncertaint y caused by the highly variable shock 
compression ratio in regions with strong hadronic particle ac celeration (e.g., iDecourchelle et alT 2 QQQI: iMiceli et al.l 
I2QQ9I: I Vink et aI]l 2 QlQl: iFerrand et al.|[2Ql5 iRessler et al.l 1201^ . We could therefore convert the measured post-shock 
plasma density to the pre-shock ambient ISM density by assuming a compression ratio of 4. This allows us to 
compare our results to the ambient ISM density inferred fro m multi-wavelength observations fe.g.. iDubner et ^120021: 
[Raymond et al .1120071: iHeng et ani2QQ7l: I Winkler et al.l[2QT^ . We caution that in our simple spectral analysis with the 
1-T plasma model, we are not able to decompose the ISM and ejecta components, which could result in an overestimate 
of the ISM density. For example, in similar X-ray measurements of the SE shell, the typ ical post-shock electron density 
is Up ^ 0.2 cm~^, cmresp onding to a pre-shock ISM density of no ^ 0.05 cm“^ (e.g., lAcero et al.l 120071: IMiceli et al.l 
1201 2 l : I Winkler et al.l 12014 ). We hereafter choose a shell thickness of 0.2 shell radius, which also produces a density of 
the NW shell roughly consistent with previous observations (ng < 1 cm“^). 

As we have adopted a simple 1 -T model, we further assume the plasma is volume filling within the X-ray emitting 
regions, i.e., / ^ 1. If the plasma is highly structured below the resolving power of our tessellated meshes, / will be 
significantly < 1. In this case, the estimated ng will be lower limit of the real value. 

From the electron density and the fitted ionization parameter ngt, we could further derive the ionization age Uon = 
netjne. In this paper, we call ngt the ionization parameter, different from some literatures where ngt is called the 
ionization age, in distinguish with tion with truly time dimension (in unit of year). As will be discussed in ^4.2.11 the 
maximum value of tj op in the SNR int erior is typically > 500 yr, consistent with the age of SN1006 (^ 10^ yr based 
on historical records: IStephensonll 20 ldl ). This also indicates that the above geometric model and our estimates of ng 
are generally reliable. 

It has been proved that the thermal X-ray emission is ejecta-dominated in some regions (e.g., the SE shell) while 
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shell thickness/shell radius 


Fig. 5.— Derived post-shock electron number density in several outer regions v.s. the assumed thickness of the thermal X-ray emitting 
shell (in unit of the outer radius of the shell). Region names are denoted in Fig. [7| 

ISM-dominated in other p arts of the SNR (e.g., the NW shell: lAcero et al.l[2QQ7l: iMiceli et alJl2QQ9Ll2Q12l:lNikQlic et alJ 
I2Q13I: I Winkler et al.l[ 2 QM ). The mixture of ISM and ejecta will not only affect the assumed abundance as discussed 
above, but also result in different shell thickness in different regions. Furthermore, the environmental density, as 
well as the filling factor of ISM and ejecta are definitely not uniform across the SNR. Therefore, we caution that the 
above estimates of Ue and Uon based on the assumption of ISM abundance, uniform shell-like geometric model, and 
volume-filling thermal plasma can be largely biased in some specific regions. 

3.3. Equivalent width map 

Equivalent width maps represent images of emission-line-to-continuum ratio, which could account for the locally 
variable underlying continuum and further reveal t he distribution of tr uly line strength throughout extended X-ray 
sources over a wide range of surface brightness (e.g.. lHwang et aT]l2QQQf ). 

3.3.1. Prominent emission lines and the eorresponding energy ranges 

In order to define the energy ranges dominated by emission lines and continuum, we present the MOS -1 spectrum 
extracted from the entire remnant in Fig. [ 6 ] (enclosed by the white solid contour in Fig. [T]). This spectrum features 
prominent emission line bumps of O, Ne, Mg, and Si, as well as weak features of S, Ar, Ca, and the K-shell blend of Fe. 
The O VII — ( lines are blended with the O VIII lines in Fig. [6l The Fe L lines often do n o t have a prom i nent b ump. 
The Fe K lines are unambiguously detected and better resolved by lYamaguchi et al.l (j2QQ8f ) ; lUchida et al.l (j2Q13f ) with 
the Suzaku observations. We define the energy ranges of the line and continuum visually and mark the emission lines 
in red in Fig. [6l The energy ranges used to construct EW maps are summarized in Table [3l We do not include the 
weak Ar and Ca lines in either the line or continuum bands. 

3.3.2. Limitations of the ahundanee map and the eontinuum estimation from linear interpolation 

In the simplest case, we obtain the continuum under the emission lines with a linear interpolation between the low 
and high energy continuums just beside the lines {Ciow and Chigh in Table [3]). The EW is calculated by subtracting 
this linearly interpolated continuum from the flux in the line band (Table [3]) and then divide the same continuum. 
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Energy (keV) 

Fig. 6.— MOS-1 spectrum extracted from the whole remnant (the white solid contour in Fig. The energy range and name of the 
emission lines are marked in red, with detailed information summarized in Table [3] “O VII high lev” means the O VII K(5 — ^ lines which 
are blended with the O VIII lines in this figure. 


TABLE 3 

Narrow bands used to created EW maps 



O VII 

O VIII 

O VII Kd -C 

Ne 

Mg 

Si XIII 

Si XIV 

SXV 

Fe K 

Line 

470-650 

650-700 

700-780 

850-1000 

1230-1450 

1660-1950 

2050-2250 

2250-2500 

6200-6700 

Clow 

300-470 

300-470 

300-470 

780-850 

1000-1230 

1450-1660 

1950-2050 

1950-2050 

4000-6200 

Chigh 

780-850 

780-850 

780-850 

1000-1230 

1450-1660 

1950-2050 

2500-2870 

2500-2870 

6700-7000 


Note. — Line, Ciow, and Chigh are the energy ranges in eV for the line, the low and high energy continuum to construct 
the EW maps, respectively. O vii K(5 — (^ is the O vii K(5, Ke, and K(^ lines which are not included in the VNEI code. These 
lines are often not clearly separated from the O viii lines in the spectra of individual tessellated regions. 


Both the abundance maps created with the spatially-resolved spectroscopy method ( ^3.1|) and the EW maps con¬ 
structed with this linear interpolation method can reveal the spatial distribution of metal-enriched gas. However, 
both methods have their own disadvantages. The accuracy of the abundance map is closely related to the goodness 
of spectral fitting. In many cases there are significant residuals at some emission lines, especially when the counting 
statistic at these lines are poor (e.g., in many cases, the Si lines are not resolved in the spectra of individual meshes). 
In addition, the tessellated meshes, constructed for spectral analysis, also have lower angular resolution compared to 
the narrow band images. 

The EW maps constructed with linear interpolations have a problem of mixing the thermal and non-thermal con¬ 
tinuums, but the emission lines are only related to the thermal emission. Therefore, it is impossible to link the EW 
calculated this way to any physical parameters such as the abundances of different elements. Eurthermore, the emission 
lines are sometimes not clearly isolated from the continuum, or the continuum shows significant curvature in the bands 
of interest. In these cases, a simple linear interpolation cannot be accurate. 

3.3.3. Continuum from spatially-resolved speetroseopy 

We herein introduce a new method to construct EW maps based on our spatially-resolved spectral analysis. We first 
calculate the flux of the continuum, i.e., both the total flux (hereafter referred as the total continuum images) and the 
flux of the thermal component only (hereafter referred as the thermal continuum images), at the energy ranges of each 
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emission lines (“Line” in Table [3|) using the fitted spectrum of each tessellated regions. We have reset the abundances 
of all the heavy elements (O, Ne, Mg, Si, S, Ar, Ca, Fe, Ni) to zero, so that the calculated fluxes roughly represent the 
underlying thermal continuum. We then use these fluxes to construct narrow band images. In order to achieve a higher 
spatial resolution, we renormalize the continuum flux within each region with the continuum images constructed with 
linear interpolations in the same band. By doing this, we have assumed this linearly-interpolated continuum images 
can roughly describe the surface brightness distribution of the total/thermal continuum within each tessellated regions. 
The total/thermal continuum images constructed this way have the same apparent spatial resolution as the original 
narrow band images. The EW maps are then constructed by subtracting the total continuum images from the original 
narrow band images, and then being divided by the thermal continuum images. 

The EW maps constructed with this spatially-resolved spectroscopy method have relatively high spatial resolution 
comparable to the original narrow band images, and are related to the thermal component only. Since we have reset 
the abundances of all the elements to zero, small residuals in spectral fitting of the line band will not affect the EW 
calculation significantly, as long as other parameters of the thermal component (e.g., the normalization, temperature, 
and ionization parameter) are well constrained. 

We caution that the presence of heavy elements will not only affect the emission lines, but also the thermal continuum. 
When the metal abundances are high, the thermal continuum could even be dominated by the heavy elements instead 
of hydrogen and helium. Therefore, the thermal-to-non-thermal continuum ratio may change significantly by setting 
the abundances of heavy elements to zero and the line images could be largely contaminated by the underlying thermal 
continuum from heavy elements. Therefore, the EW map constructed this way should contain the contributions from 
both emission lines and the thermal continuum produced by heavy elements, and should be adopted as an upper limit 
of the true EW. 


4. RESULTS AND DISCUSSIONS 
4.1. Speetra from larger regions 

In order to check the reliability of parameter maps presented in the following sections, we extract spectra from larger 
regions (than the tessellated meshes) enclosing some interesting features. These regions are denoted in Eig. Oand the 
spectra extracted from them are presented in Eig. [51 We calculate th e ave rage values of some parameters of the thermal 
component based on the parameter images as will be presented in ^4.21 in order to roughly characterize the average 
properties of each regions. Eurthermore, we also fit the sp ectra extracted from these larger regions using the same 
model as adopted for the much smaller tessellated meshes ( ^3.2.ip . in order to compare with the average parameters. 
We caution that there is generally a poor fitting of these high-counting-statistic spectra with this simple 1-T thermal 
plasma model, and the fitted parameters are only valid for a general comparison of the trend among different regions. 
The averaged parameters and the parameters obtained from spectral analysis are both summarized in Table jH 

As presented in Table IH the fitted parameters are typically lower than the averaged parameters, especially for 
logTT-et, and the O, Ne, Mg, Si abundances. This is primarily caused by the extreme values of the parameters in some 
meshes with inadequately fitted spectra. As shown on the parameter maps (e.g., Eig. [9]), these extreme values are 
just in a few isolated meshes, but could significantly bias the average value of certain parameters in a large region. 
Eor example, the “NW Shell” and “SNR Interior 01, 02” regions have several white pixels on the Uet map (extremely 
large Uet values of > 10^^ cm“^s; Eig.jHb), so their averaged log net show the largest differences from the fitted values. 
In conclusion, except for some extreme values (or “bad pixels”), the parameter maps constructed with our spatially 
resolved spectroscopy method are generally consistent with the spectral analysis results of larger regions. 

We use three different ways to describe the spatial distribution of emission lines from different ions: (1) the abundance 
map constructed with spatially resolved spectroscopy (hereafter 2-D Spec method); (2) the EW map constructed with 
either linear interpolation of the continuum (hereafter linear EW method) or (3) the continuum measured from the 
thermal component of the spatially resolved spectral modeling (hereafter 2-D Spec EW method). In principle, the 2-D 
Spec method is the most reliable in physics but has the lowest spatial resolution. In contrast, the linear EW method 
has the highest spatial resolution but cannot be directly linked to physical parameters because of the contribution 
from the non-thermal component and the curvature of the soft X-ray continuum ( ^3.3.2p . The 2-D Spec EW method 
represents a compromise of these two methods, with the continuum determined from the 2-D Spec method, while the 
line images constructed in a similar way as for the linear EW method ( ^3.3.3|) . As the O abundance may be affected 
by O VII, O VIII, and O VII KS — ( lines, we present in Table |4] the EWs (from both linear EW and 2-D Spec EW 
methods) of these three lines together with the O abundance. The relative strength of these lines also indicates the 
thermal and ionization states of the plasma. 

4.2. Parameter maps from spatially-resolved speetroseopy 

We present the parameter maps constructed with the spatially resolved spectroscopy method in Eig. |9l In the 
following sections, we describe interesting structures on each map and discuss their physical implications. 

4.2.1. Thermal and ionization states 

The thermal and ionization states (characterized by kT and n^t) trace the shock heating and ionization history of 
the S NR. In youn g SNRs, the forward and reverse shocked gas is often far from temperature and ionization equilibrium 
fe.g.. lSlan 3 l20l4 . X-ray measurements in the literatures reveal a large diverse of kT and n^t in SNI006 (with either 
I-T model or more complicated models), with kT ranging from < 0.5 keV to > 4 keV, and Uet typically in the 
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Fig. 7.— Regions used to extract the spectra presented in Fig. [8] and summarized in Tabled with region names denoted beside. The 
background tricolor images are the same as those presented in Fig. [T] 
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Fig. 8.— Top row (panels a-d): the MOS-1 spectra of the 12 regions (Fig. [7|) with the region names denoted on the top right corner of 
each panel. We present the spectra in four different panels to clarify. Bottom row (panels e-h): similar as the top row, but for the PN 
spectra. MOS-2 spectra are similar as MOS-1, so are not presented here. For the convenience of comparison between the MOS and PN 
spectra, each spectra are divided by their effective area (in unit of counts s“^ keV“^ cm“^ instead of counts s“^ keV“^). The x and y 
axes of each panel are in the same range. 


rang:e of a fewxlO^ cm~^s fe.g.. iKovama et ^119951: IVinket apl2QQQl: iDver et alll^OOll: lAllen et alll2QQll: IVink et ^ 
I2QQ3I: iLong et aI]l2QQ3l: lAcero et al.ll2QQ7l: lYamaguchi et al"]l2QQ8l: iMiceli et alJl2Q12l: lUchida et alJl2Ql^ . Although the 
presence of non-thermal emission in the NE and SW synchrotron rims may lead to large systematic uncertainties in the 





























13 


TABLE 4 

Average value of parameters for individual regions 


Region 

Method 

kT 

keV 

rie 

cm“^ 

log net 
log(cm“^s) 

O VII KS-C 
EW 

O VII 
EW 

O VIII 
EW 

O 

solar 

Ne 

solar 

Mg 

solar 

Si 

solar 

Fe 

solar 

NW Shell 

Average 

Eit 

2.25 

1.58 

0.52 

0.39 

9.56 

9.32 

0.18, 23.3 

0.52, 12.0 

0.33, 51.9 

1.13 

0.92 

0.85 

0.64 

2.20 

1.05 

11.06 

4.84 

0.23 

0.05 

SNR Interior 01 

Average 

Fit 

2.10 

2.22 

0.26 

0.24 

9.50 

9.37 

0.15, 33.1 

0.57, 15.3 

0.33, 61.9 

1.77 

0.90 

0.49 

0.27 

3.81 

1.14 

32.38 

7.85 

0.18 

0.39 

SNR Interior 02 

Average 

Fit 

1.63 

1.56 

0.33 

0.31 

9.51 

9.35 

0.21, 38.3 

0.64, 17.8 

0.41, 74.9 

2.07 

1.05 

0.98 

0.44 

5.11 

1.69 

30.71 

12.36 

0.52 

0.83 

SNR Interior 03 

Average 

Fit 

2.47 

4.59 

0.34 

0.28 

9.40 

9.35 

0.17, 23.8 

0.60, 13.2 

0.36, 51.5 

1.44 

1.23 

0.50 

0.32 

3.42 

1.46 

6.95 

5.31 

0.25 

0.05 

SNR Interior 04 

Average 

Fit 

2.87 

2.30 

0.37 

0.35 

9.48 

9.40 

0.19, 27.0 

0.60, 10.6 

0.44, 52.3 

1.26 

1.01 

0.36 

0.28 

1.97 

1.53 

9.83 

7.36 

0.40 

0.36 

SNR Interior 05 

Average 

Fit 

1.30 

1.12 

0.41 

0.36 

9.65 

9.60 

0.21, 34.8 

0.43, 10.8 

0.37, 59.3 

1.53 

0.80 

0.62 

0.35 

3.16 

1.26 

20.23 

9.31 

0.18 

0.28 

Dark Belt 

Average 

Fit 

2.06 

2.54 

0.43 

0.32 

9.49 

9.41 

0.17, 25.1 

0.46, 8.99 

0.33, 44.7 

1.06 

0.87 

0.39 

0.32 

1.10 

0.94 

15.33 

9.30 

0.73 

0.68 

Interior Shell 01 

Average 

Fit 

2.41 

2.20 

0.54 

0.43 

9.58 

9.53 

0.22, 35.9 

0.40, 8.43 

0.37, 50.5 

1.06 

0.97 

0.42 

0.38 

0.86 

0.82 

9.61 

8.86 

0.99 

0.96 

Interior Shell 02 

Average 

Fit 

2.36 

1.61 

0.53 

0.51 

9.59 

9.54 

0.20, 31.1 

0.37, 8.85 

0.37, 53.0 

1.23 

0.77 

0.35 

0.26 

1.62 

1.01 

13.41 

8.67 

0.34 

0.44 

O hole 

Average 

Fit 

1.70 

1.65 

0.51 

0.55 

9.61 

9.61 

0.16, 31.1 

0.28, 5.74 

0.27, 38.7 

0.81 

0.43 

0.30 

0.21 

0.84 

0.47 

13.93 

5.48 

1.04 

0.63 

SE Shell 01 

Average 

Fit 

1.89 

2.32 

0.53 

0.42 

9.53 

9.45 

0.17, 28.3 

0.40, 8.38 

0.37, 47.5 

1.10 

0.99 

0.27 

0.24 

1.44 

1.32 

16.56 

12.90 

0.69 

0.76 

SE Shell 02 

Average 

Fit 

1.39 

1.07 

0.46 

0.41 

9.52 

9.58 

0.16, 27.5 

0.51, 10.0 

0.39, 53.2 

1.33 

1.23 

0.34 

0.31 

2.65 

1.82 

22.06 

17.61 

0.28 

0.18 

SE Shell 03 

Average 

Fit 

1.65 

1.66 

0.32 

0.30 

9.54 

9.48 

0.22, 27.4 

0.64, 9.63 

0.60, 55.0 

1.45 

1.15 

0.28 

0.22 

3.21 

2.24 

15.85 

11.52 

0.10 

0.05 


Note. — Average parameters of large regions enclosing some interesting features as denoted in Fig.[3 For each region, the average parameters 
are calculated in two ways: a direct average based on the parameter images (“Average”) and the parameters obtained by fitting the MOS-l+MOS- 
2+PN spectra extracted from each region (e.g., Fig.[^ using the model described in asunK “Fit”). For the former method, /cT, log net, and Ue 
are calculated from Fig. it, b, and d. O vii, O viii, and O vii K5 — 4 EWs are calculated from the linear EW maps presented in Fig. fT2k-c 
(former numbers) and the 2-D Spec EW maps presented in Fig. I12t i-f (latter numbers). O, Ne, Mg, Si, and Fe abundances are calculated from 
the abundance maps in Fig. I12h . Fig. 113b . Fig. 114b . Fig. 115b . and Fig. 117b . 


derived kT and net, which may partly account for the spread in the reported kT and Uet measurements, for the interior 
of the SNR the range of measured kT and Uet likely reflect intrinsic temperature and ionization state variations. 

In the present paper, we adopt a simple 1-T model to characterize the thermal plasma emission ( ^3.2.ip . kT and 
Uet are thus highly affected by the most prominent O lines, as well as the Ne, Mg, and Fe L lines, but not significantly 
affected by the poorly resolved Si and S lines. We present the kT and Uet maps in Fig.[9^,b, and caution that kT and 
Uet in the NE and SW non-thermal limbs are not well constrained and could be significantly overestimated in some 
regions. 

In the regions dominated by thermal emission, the thermal and ionization states of the plasma show clear spatial 
variations. kT typically varies in the range of ^ I — 4 keV, while Uet typically varies in the range of 10^“^^ cm“^s. 
These ranges are generally consistent with previous measurements. The outer shells of the SNR (the NW and SE 
shells) appear to have low kT and Uet, indicating that the gas there are newly shocked and are far from thermal and 
ionization equilibrium. kT of the SNR interior is relatively smoothly distributed, with most of the regions having 
kT ^1.3 keV. Some higher temperature structures are superimposed, such as “SNR Interior 03 and 04”, “Interior 
Shell 01 and 02”, etc. These features are typically locally high-density structures. In contrast, the Uet map shows clear 
gradient in spatial distribution, with the highest Uet appears close to the center of the SNR (“O hole”), and declines 
outward. The “Dark Belt” has clearly lower Uet than the surrounding regions. 

We characterize the distribution of the thermal and ionization states with the original and emission-measure-weighted 
(EM-weighted) probability distribution functions (PDEs) (Eig. [10]). The EM-weighted PDEs of both kT and Uet can 
be characterized with a 5-Gaussian model. Interestingly, the three peaks at 0.22 keV, 1.32 keV (or the broad Gaussian 
centered at 1.59 keV), and 4.35 keV in the kT PDE are roughly consistent with the three components from the spectral 
analysis of lUchida et al.l (j2013f ). where they attribute the lowest temperature component to shocked ISM, while the 
other two components to shocked ejecta. Nevertheless, as shown in Eig.flOb.c. three of the five Gaussians are relatively 
narrow and have low EM (the kT = 0.22, 4.35, and 10.6 keV components); they are either in the non-thermal emission 
dominated regions so have poorly constrained kT and Uet or are in some small isolated regions in the SNR interior. 
We therefore only use the two most prominent Gaussians to characterize the average thermal and ionization states of 
SN1006. We assume the kT = 1.32 ± 0.37 keV component and the log net/(cm“^s) = 9.48 ± 0.14 component are from 
the same regions, while the kT = 1.59 ± 1.61 keV component and the lognet/(cm“^s) = 9.84 ± 0.47 component are 
from the other regions. 

We further compare our measured kT and Uet of all the tessellated regions to those obtained from the archival 
measurements of SN1006 (Eig. [TT]) . We also plot the two most prominent components in the PDEs of kT and Uet to 
represent the average thermal and ionization states of SN1006. The archival measurements are obtained with different 
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Fig. 9. — Parameter images from spectral analysis of individual tessellated regions as shown in Fig. [21 with color bar on right and 
parameter name on the top right corner of each panel, (a) Temperature of the VNEI component in unit of keV (kT); (b) Ionization 
parameter of the VNEI component in unit of cm“^ s (net); (c) Ionization time in unit of year (tion)] (d) Electron number density in unit 
of cm“^ (^e); (e) Radio spectral index of the SRCUT component (o; ); (f ) Cutoff frequency of the SRCUT component (i^cutoff)- White 
regions overlaid on each panel are denoted in Fig. [7| and described in fl4.1l 
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Fig. 10. — Probability distribution functions (PDFs) of the temperature (kT) and ionization parameter (net) of the VNEI component of 
all the 3596 tessellated regions. Panel (a) is the PDF of kT with the y-axis denoting the number of regions, while panels (b) and (c) are 
the emission-measure-weighted PDFs with the y-axis denoting the total emission measure. The solid curve in panels (b) and (c) are the 
best-fit multi-Gaussian model, with the dashed curves mark the model components. The central positions of each Gaussian component are 
also denoted. The errors are at 1-a level, obtained from the width of the Gaussian component. 


data, from different regions, and/or with different models, so span a large range of kT and Uet. Nevertheless, most of 
these measurements are roughly consistent with our measurements. Therefore, the large diversity in the measured kT 
and Uet from previous works may simply because they are obtained in different regions consisting of plasma at various 
thermal and ionization states. 

Combining the Ue and Uet maps, we construct the ionization age {tion) map with truly time dimension (in unit of 
year; Fig. |9t). Except for the bright rim surrounding the SNR, which is artificial due to the low flux density of the 
surrounding regions, the whole SNR shell appears to have a low and smooth ionization age of tion ^150 year. In 
contrast, all the regions in the SNR interior have tion >150 year, with the highest tion ^ 500 year, consistent with 
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Fig. 11.— kT v.s. net for all the 3596 tesselllated regions (black dots). The red filled circles (“PDF”) are the two most prominent 
Gaussian components in Fig. llOb .c: the strongest component with kT = 1.32 ± 0.37 keV and lognet/(cm“^s) = 9.48 ± 0.14 and the 
broadest component with kT = 1.59 ± 1.61 keV and log net/(cm“^s) = 9.84 ± 0.47. Other colored symbols are archival measurements of 
SN1006 with the references denoted in the lower left corner. 
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the age of SN1006 of ~ 10^ year. The tion map suggests that the plasma in the SNR shell is probably newly shocked, 
while the SNR interior was shocked several hundred years before after the SN explosion and is on the approach of 
ionization equilibrium. 


4.2.2. Electron density 


Many previous studies indicate that SN1006 locates in a low-density environment. The measured density of the 
surrounding ISM is typically in the range of 0.04 ^ no ^ 0-4 c m~^ and the NW pa rt is evidenced to have th e highest 
densi ty. Different estimates of the ISM density are based on H I (iDnbner et~^l 2002 D. IR (IWinkler et al"]l2013[ ). Balmer 
lines (jRavmond et al.ll^OOTHNikolic et al.ll2013D . and UV (iLaming et al.lll996l:lKorreck et 2004) observations, as well 
as th e thermal component of the X-ray emission (j Winkler fc Longill997l : lLong et al.ll2003l: lAcero et al.ll2007l iMiceli et al.l 
l2012f ). This low environmental density also explains SNlOOO’s faintness in TeV 7 -ray emission, which could partly 
come from the proton-proton interaction with tt^ production and subsequent decay (jAharonian et al.l[2005l: lAcero et~al] 

[2oTol . [2oTa . 


For the first time, we obtain the map of electron number density of the thermal plasma (ue) over the entire SNR. 
As shown in Fig. [9]i, the overall Ue distribution resembles the total surface bright ness distribution, ex cept for the two 
non-thermal limbs. The Ue image also resembles the “thermal” image presented in IMiceli et al.l (j2QQ9l ). consistent with 
its origin from the thermal emission. The “NW shell” has significantly higher Ue than the “SNR Interior”. Assuming 
no CR acceleration in this part so a compression ratio of 4, the ambient ISM density surrounding the “NW shell” 
should be no ^ 0.15 cm“^, consistent with previous multi-wavelength estimates. There is no sharp increase of Ue in the 
SE part. Instead, the “SE Shell” appears much thicker than the “NW Shell”, and the SE half of the SNR is separated 
by the low density “Dark Belt” into two shells (the “SE Shell” and the “Interior Shell”), both with comparable Ue 
as the “NW Shell”. The “Dark Belt” has low Ue and Uet compared to the surrounding regions (“Interior Shell”, “SE 
Shell 01”), but has tion comparable to or even higher than these regions. The low ionization states of the ions in the 
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“Dark Belt” are further indicated by its relatively low O VIIl/0 VII ratio (Fig.[8t; Table S]) and low centroid energy of 
the Si XIII lines (Fig. [St, g). The outermost part of the “SE Shell” (e.g., “SE Shell 03”) decreases in ng. The northern 
part of the “SNR Interior” has a clearly lower density of Ue ^ 0.3 cm“^, but “SNR Interior 03 and 04” may form 
a shell-like structure behind the SW non-thermal rim, apparently extending the “NW Shell”. In higher resolution 
intensity images (e.g., Eig. o, many parts of the SNR interior, such as the “SE Shell” and the “Interior Shell”, are 
resolved into clumpy structures, in contrast to the significantly filamentary structure in the “NW Shell”. 

Based on the Ue map, we could further estimate the total mass of the plasma contained in the 3596 tessellated regions, 
again assuming a solar abundance. The estim ated m ass of the shocked X-ray emitting plasma is ^ 14 . 5/2 M©, where 
/ is the volume filling factor. As discussed in 93.2.21 the swept up ambient ISM mass is ^ 5 M©. Adding the mass of 
the shocked ejecta, which is quite uncertain but contribute only a small fraction to the mass budget, we could roughly 
constrain the volume filling factor to be / ^ 0.1 under the adopted geometric model f 93.2.2p . 

4.2.3. Non-thermal component 

The shape of t he non-thermal X-ray spectra is an indicator of many parameters of particle acceleration at the 
SNR shocks (e.g.. [Alien et al.ll2QQlL \200A [Rothenflug et al.[r2QQ4 [Miceli et~ffl[2QQ9L [20131 ). We will further discuss the 
modeling of the non-thermal spectrum of SN1006 in companion papers, including our newly approved NuSTAR Cycle 1 
observations of the non-thermal rims of SN1006 (400 ks; PI: Li Jiang-Tao). We herein present the preliminary results 
on the non-thermal X-ray emission from the simplest spectral decomposition with the “VNEI-hSRCUT” model. As 
descr ibed in 93.2. 1[ we fi x the normalization of the SRCUT component using the 1.4 GHz flux-accurate image of SN1006 
from [Dyer et al.[ (|2009[ ). Therefore, the only two free parameters of the SRCUT component are the radio-to-X-ray 
photon index a and the cutoff frequency I'cutoff^ which are presented in Eig. [9fe,f. 

a of the regions dominated by non-thermal emission is quite smooth, roughly in the range of 0.4-0.5. This is 
consi s tent with the adopti on of a constant a ^ 0.45 — 0.6 in some previous works (e.g., [Dver et ^[2001[: [Acero et al.[ 
[2007[: [Katsuda et al.[[2010[ ). However, it is also possible that there is a more complicated shape o f the non-thermal 
X-ray emission an d a is varying not only at different locations, but also at different energies (e.g., [Allen et al.[ [20081: 
[Miceli et al.[[2013[ ). The current energy coverage (with Chandra or XMM-Newton) is too narrow in hard X-ray and we 
expect a better modeling of the electron synchrotron emission with our new NuSTAR data. 

The cutoff frequency nn y.tn f f shows significant s p atial variations with in the non-thermal emission dominated regions, 
as previously revealed by lR.othenfiug et al.l (j2004f ) ; [Miceli et al.[ (j2009[ ) . The bright non-thermal filaments are resolved 
on the ^cutoff map and have higher ^cutoff of > 5 x 10^^ Hz, while the inter-filament regions typically have ^cutoff ^ 
(2 — 5) X 10^^ Hz. ^cutoff continuously falls toward the SNR interior, and is < 5 x 10^^ Hz in the regions dominated 
by thermal emission. 


4.3. Spatial distribution of metals 

Altho ugh the presence of emission lines from metal ions in SN1006 has been suggested in some early X-ray obser- 
vatio ns ([Becker et al.[ [198Q[: [Leahv et al.[[l99lh . these lines are not reliably detected until the publication of ASCA 
data (jKovama et al.[[l995[ ). Later on, with BeppoSAX^ Chandra^ XMM-Newton^ and Suzaku observations, the emis¬ 
sion lines are better resolved and the spatial distributions of metal ions have been extensively st udied, either with 
spectral analysis of different regions or the narrow band uiapping of different emissi on lines (e.g., [Vink et ^ [2QQQL 
[2QQ3[: [Long et ^ [2QQ3[: [Acero et~aLl[2QQ7[: [Yamaguchi et al.[l2QQ8[ : lUchida et al'][2Q13[ ). We herein study the spatial 
distribution of different elements in SN1006, using both the abundance map and the EW maps constructed with the 
linear EW or 2-D Spec EW method. 


4.3.1. O 

There are usually two prominent oxygen line bumps presented in the soft X-ray spectrum of SN1006: the O VII 
K-shell transitions at ^ 0.56 — 0.58 keV and the O VIII K/3 and Ky lines at ^ 0.65 keV, and sometimes a third line 
at ^ 0.7 — 0.75 keV representing the O VII KS — ( transitions (Eig. [6]). In Eig. [121 we present the O VII, O VIII, and 
O VII KS — ( EW maps constructed with both the linear EW method and the 2-D Spec EW method, the 2-D Spec 
EW ratios of these lines, as well as the O abundance map. 

As shown in Eig. [121 fhe intensities of O VII and O VIII gradually decrease toward the center, with the lowest 
intensity at the “O Hole”. This spatial variation indicates a gradual decrease of the O abundance toward the center. 
The O abundance across the whole SNR is close to solar value, except for the non-thermal arcs where the properties 
of the thermal component are not well constrained. Therefore, expect for some fine structures, such as the “SE Shell 
02 and 03” (also refer to Table [4]), we generally do not find significant evidence of O enrichment in SN1006. 

The O VII and O VIII distributions show some noticeable differences. As shown on the O VIIl/0 VII map (Eig. [121), 
in the NW edge, the O VII lines appear to be relatively stronger, while the O VIII lines are stronger in the opposite 
side (the SE half of the SNR). The O VIIl/0 VII map resembles the Uet map (Eig.[9l)), indicating that Uet is primarily 
determined by the most prominent O lines. As revealed by the O VIIl/0 VII and Uet maps, the plasma in the SE half 
of SN1006 is closer to ionization equilibrium. However, this higher Uet is most likely a result of the higher density in 
this region, as the ionization time {tion) at the SE and NW edges seem to be comparable (Eig. [9t). 

The only feature with significantly subsolar O abundance is the “O Hole” (Table [4]). It is interesting that this 
region, in addition to the “Interior Shell 01” close to it, also has significantly enhanced O VII KJ — C/0 YII and 
O VII KS — C/0 VIII ratios (Eig. [12^, h). The O VII KJ — C/0 VII map is also very similar in morphology as the Uet 
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Fig. 12.— (a-c) O VII, O VIII, and O VII K(5 — C EW maps constructed with the linear EW method. Background pixels are set to 0, so 
pixels with color bluer than the background color have negative values, (d-f) Similar as panels (a-c), but constructed with the 2-D Spec EW 
method, (g-i) The EW ratio of O VII K(5 — C/O VII, O VII K(5 — C/O VIII, and O VIIl/O VII constructed with the 2-D Spec EW method, 
(j) Oxygen abundance map obtained by fitting the spectra in each tess ellat ed regions, with the color bar in unit of solar abundance. White 
regions overlaid on each panel are denoted in Fig. [7| and described in su 


map, indicating that the regions with strong O VII K(5 — ( emission apparently have high ionization state. However, 
this may instead be an artifact because it is easy to mix O VII K(5 — ( and O VIII lines. Alternatively, the O VII K(5 — ( 
lines may also be highly contaminated by the low ionization transitions of Fe L lines. We will further discuss this 































18 


possibility in 114.3.51 


4.3.2. Ne 

In case of no prominent iron L-shell lines, the significant emission line bump at ^ 0.9 keV mainly consists of Ne IX 
K-shell transitions (Fig. [6l Table [3]). Some regions with strong Ne emissions also show weak Ne X K-shell emission 
lines (at ^ 1.0 keV; e.g., the “NW Shell”; Fig. [8^), which are much less significant in the spectra of other regions. 

As shown in Fig. [131 fhe Ne abundance map (panel c and d) shows strikingly similar features as the Ne EW maps 
(panels a and b). Some features with enhanced Ne emission are clearly resolved, such as the “NW Shell”, the weak 
enhancement roughly at the “Dark Belt” (or a little behind, i.e., toward NW) and some parts of the “SNR Interior”. 
In contrast, the “SE Shell” has extremely low Ne emission. The coincidence of these features in different panels of 
Eig. Uni strongly suggest the enhancement at the NW shell is related to the SNR itself, instead of just an artificial 
background feature. 
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Fig. 13. — The Ne IX EW maps constructed with (a) the linear EW method and (b) the 2-D Spec EW method, (c) The Ne abundance 
map obtained by fitting the spectra in each tessellated regions, (d) Ne/O abundance ratio map. 


Except for the “NW Shell” (Ne/O > 0.7; Table |4]), the Ne/O ratio is significantly below solar across the entire 
SNR (Fig. [T3lib The highest Ne/O ratio inside the SNR is < 0.5 (at the “Dark Belt” and “O Hole”) and the typical 
value of the Ne/O ratio is ^0.2-0.3. The NW high-Ne filament also has high electron density, as well as strong soft 
X-ray and Ha emissions, which are often signatures of the shocked ISM (e.g.. IWinkler et al]l2014l ). Therefore, Ne in 
SN1 Q06 is most likely fro m the ISM instead of the ejecta, consistent with the small amount of Ne yielded in Type la 
SN (|Iwamoto et al.l[l9^ . If the high-Ne regions represent shocked ISM, it is very likely that a signific ant fra ction of 
the O line emissions are from the SN ejecta, although the O abundance is mostly solar across the SNR ( ^4.3.ip . as the 
O lines are always weak in high-Ne regions. 


4.3.3. Mg 

In case of no prominent iron L-shell lines, the significant line bump at ^ 1.3 keV mainly consists of Mg XI K-shell 
transitions and possibly some weak Ne X K-shell lines (Fig. |6l Table [3]). Similarly, the Mg abundance map also shows 
coincident features to the Mg EW maps constructed with either the linear EW or the 2-D Spec EW method (Fig. [T4|) . 
such as the low intensity at the “O Hole”, “Interior Shell 01”, “Dark Belt” and the relatively high i ntensity in the 
“SNR Interior” and “SE Shell”. Type la SNe are not expected to yield a lot of Mg (|Iwamoto et al.l[l999f ). so Mg should 
be mostly ISM in origin. However, the spatial distributions of Mg and Ne are quite different. For example, compared 
to the Ne lines, the Mg line strength is pretty low at the “NW Shell” but high at the “SE Shell” (Fig. 0 Table |4]). 
These differences in Ne and Mg distributions indicate that either SN1006 is Mg enriched due to some unknown reasons 
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or it is more likely that the ISM conditions (density, ionization state, etc.) are significantly variable around SN1006, 
which produces variable Mg/Ne emission line ratios even though their abundance ratio is not significantly variable. In 
most of the regions, the Mg/0 ratio is > I (Fig. [Mhb However, this is very likely to be an artificial effect as we have 
adopted a I-T model which often fit the most prominent O lines well so may have overpredicted the abundances and 
EWs of emission lines at higher energy. 
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Fig. 14.— Similar as Fig. 1131 but for Mg XI lines and the Mg abundance. 
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4.3.4. Si and S 

Because of the low counting statistics. Si lines are often weak in the spectra of individual tessellated regions, and 
S lines are even weaker. We therefore link the Si and S abundances in spectral fitting, by assuming that they are 
produced in related processes in stellar nuclear synthesis (e.g., oxygen burning) so are co-spatial. There are two 
prominent Si line bumps in the spectra of SN1006 (Fig.[ 6 |). The Si XIII bump centers at ^ 1.8 keV and mainly consists 
of K-shell transitions of helium-like silicon ions, although some higher level transitions and Mg XII K-shell lines may 
also contribute. There is also an emission line feature centered at ~ 2.1 keV, and is most likely blueshifted Si XIV 
and Ky lines with a rest frame energy of ^ 2.0 keV. Another weak bump centered at 2. 4- 2. 5 keV is most likely S XV 
K-shell transitions, although some weak Si XIV or even Si XIII K-shell transitions may contribute. 

The Si XIII and Si XIV EW maps and the Si abundance map are presented in Eig.[T5l The Si XIV lines are very weak 
and have no detectable high signal-to-noise ratio features on the images, although some enhancements in the SNR 
interior may still be possible. On the other hand, the Si XIII lines are quite prominent and determine the measurement 
of Si abundance. The Si XIII line bump is the strongest in the “SE Shell” and shows clearly asymmetric distribution. 
The outermost shells (“NW Shell” and “SE Shell 03”), as well as the regions slightly NW to the “Dark Belt”, appear 
to be relatively weak in Si XIII emission. The Si/0 ratio is above solar (Eig. fTbt). Similar as the Mg/0 ratio, it may 
also be an artificial effect caused by the 1-T thermal plasma model. 

“SNR Interior 01-05” represent regions selected from the low surface brightness interior between the NW shell and 
the relatively bright SE half. The five regions have the same surface area, so their spectra shown in Eig. [H^,b and e,f 
are also indicative of the intrinsic intensity of each region. It is clear that “SNR Interior 01, 02 , and 05” are fainter 
than the other two, which apparently coincide with two shell-like features projected just behind the west non-thermal 
shell. “SNR Interior 01” has the lowest electron density, consistent with its low surface brightness. Compared to “SNR 
Interior 01, 02 , and 05”, the Si abundances (or the Si/0 ratio) of “SNR Interior 03 and 04” are lower (Table |4]), as also 
indicated by their weaker Si lines (Eig. [H]3,f)- It is therefore likely that “SNR Interior 03 and 04” are largely comprised 
of shocked ISM if the Si lines are mostly from the SN ejecta. 

In order to study the possibly different distributions of Si and S, we also present the S XV EW maps constructed 
with the linear EW method and the 2-D Spec EW method in Eig. [161 The S XV EW map constructed with the 2-D 
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Fig. 15.— Panels (a-c) and (f) are similar as Fig. 1131 but for Si XIII lines and the Si abundance. Panels (d, e): the Si XIV EW maps 
constructed with the linear EW method and the 2-D Spec EW method. 
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Fig. 16.— The S XV EW maps constructed with (a) the linear EW method and (b) the 2-D Spec EW method. 


Spec EW method is quite noisy, probably because of the uncertainty in modeling the spectra at low signal-to-noise 
ratio. The S XV EW map constructed with the linear EW method better traces the S distribution, which is quite 
similar as the Si XIII lines. The smoother appearance is apparently a result of the lower signal-to-noise ratio of S lines. 

4.3.5. Fe 

Type la SNRs are oft e n thought to be the major source of Fe in t he hot g as in and around galaxies (e.g., 
iHnmphrev fc Bnot3l2QQ6l : iJi et al.lI2QQ9I: iLi et al.l 120091 : iLi fc Wan3l2013l : IlTI 120151 ). However, as a typical Type la 
SNR, SN1006 is well known to exhibit a “missing iron” problem, i.e., has no significant Fe emission lines in most 
parts of t he SNR fe.g..lVink et al.ll2003r). Neverthel e ss, th e existence of Fe is still revealed in some different ways. For 
example. lYamagnchi et aP (|2008f ) and lUchida et al.l (|2013f ) have recently reported the detection of Fe K emission lines 
(at a very low ionization state of riet < 10^ cm“^s and center at ^ 6.4 keV) with deep Suzaku observations. On the 
other hand, optical and UV observations of brigh t background sources toward the direction of SN1QQ6 have revealed 
some Fe absorption hues at low ionization states (jFesen et al.lll98^ IWn et aH 119931: iBlair et allll996l: iHamilton et al.l 
Il997l : IWinkler et al.ll2QQ5] ). Therefore, it is most likely that the Fe-rich layers of SN1006 has not yet been completely 
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shock heated by the reverse shock. In fact, these UV/optic al spectra reveal that the reverse shock is penetrating still 
in a Si-rich zone (j Winkler et alJ[2QQ^ iHamilton et alJ[2QQ7[ ). 



Fig. 17.— (a) The Fe K line EW map constructed with the linear EW method, (b) Fe abundance map. (c) Fe/O abundance ratio map. 
We caution that the Fe abundance is mainly determined with the widely distributed Fe L lines at ~ 1 keV, different from the Fe K line 
shown in panel (a). The Fe L lines typical do not appear as a prominent bump. We therefore do not construct EW maps for them. 


Spatial distribution of Fe is shown in Fig. [T^ The Fe K-shell lines are very weak and only marginally detectable in 
our XMM-Newton spectra extracted from the entire remnant (Fig. [6]). In addition, the thermal continuum is mainly 
determined in soft X-ray (< 2 keV) and is thus poorly constrained in the Fe K line band 6.2 — 6.7 keV). Therefore, 
we only present the Fe K-line EW maps constructed with the linear EW method in Eig. flTh . together with the Ee 
abundance maps (Eig. [TTb.c) mainly determined with the Ee L-shell lines at ^ 1 keV. It is obvious that the strongest 
Ee emission appears in the SE quarter, consiste nt with the narrow band image and the sp atially resolved spectral 
analysis conducted with the Suzaku observations (|Yamaguchi et all|2 0081: llJchida et al.ll2Q13f) . 

It is interesting that the Ee abundance distribution (Eig. flTb.c) is quite similar as the O VII KS — C/0 TW 
ratio (Eig. [T2]i). The peak around the “O Hole” is so sharp that the O VII — C/0 VIII ratio is at least twice 
lower and there is negligible Ee emission in other regions. The presence of the 6.4 keV Ee K shell lines, as well as 
the relatively strong O VII KS — C lines, indicate that the ionization parameter should be low around the “O Hole”. 
However, the measured Uet peaks at this region (Eig.lHb) apparently lead to opposite conclusio ns. Therefore, there are 
at least two compo nents of the plasma toward the “O Hole”, one with low Ugt of < 10^ cm“^s (|Yamaguchi et al.ll2QQ8l: 
lUchida et al.ll2Ql^ accounting for the Ee K and O VII KJ —C emission lines, the other with high Uet of ^ 5 x 10^ cm“^s 
accounting for most of the soft X- ray emission includin g the O VII and O VIII lines. The high-riet component is most 
likely from the shocked ISM (e.g.. lEerrand et aIll2Q12f ). which is also consistent with the roughly solar O abundance. 
On the other hand, the low-Uet component is probably the reserve shocked ejecta, which also explains the high Ee 
abundance. In the Type la SNR SN1006, only the central region shows significant Ee emission (from both K-shell and 
L-shell transitions), which means the reverse shock has just reached the Ee-rich shell of the ejecta after > 10^ yr since 
t he SN explosio n. _ 

lUchida et al.l (|2Q13[ ) suggested another possibility to explain the spectral fitting residual at 0.7-0.75 keV. There could 
be relatively strong Ee L 3d ^ 2p transitions of extremely low-ionization Ee ions (the most dominant charge state of 
Fe is l ess than Fe^^+), in addition to the O VII Kd — C lines. This possibility is also suggested bv IWarren fc HnghesI 
(|2004D for the Type la SNR E0509-67.5. Such contributions from newly shocked low Fe ions are also evidenced by the 
elevated emission at ^0.8-0.9 keV of “SE Shell 01” (Fig. [Hli,h), which has similar properties as “SE Shell 02 and 03” 
but has significantly higher Fe abundance (Table 01 Fig. [TTb.cb The Fe L 3d ^ 2p transitions of low Fe ions naturally 
explains the spatial coincidence of the 0.7-0.75 keV residual and the Fe abundance. 

Many geometric models of SN1006 have been proposed in the literature, in order to explain t he asymmetric distri¬ 
bution of multi-wav elength features, primarily based on UV/optical absor ption line studies (e .g. jHarndton et al.lll997l: 

I Winkler et al.l [200^ or X-ray mapping of prominent emission lines (e.g., lUchida et al.ll201^ . In the above sections, 
we for the first time present maps of many physical parameters of SN1006, such as /cT, rie, Uet, and Uon- We also 
present abundance and EW maps of many heavy elements. Most of these parameter maps show clear asymmetric 
distributions, which could be attributed to the asymmetric distributions of either the ambient ISM or the SN ejecta. 
For example, the Ne emission lines are very stron g in the NW shell, consistent with an e nhanced ISM density in this 
region as revealed in H I and Ka observations (e.g. dDnbner et al~]l2QQ2l: iNikolic et al.ll2Ql^ also see the electron density 
map shown in Fig. [9li). However, the very different distributions of O, Ne, and Mg emission lines, which are usually 
thought to be good tracers of shocked ISM, indicate that the asymmetric distribution of ISM cannot explain all the 
characteristics of the spatial distribution of metals. The very asymmetric distribution of Si indicates that the reverse 
shock encounters the ejecta earlier in the SE half of the SNR. The detection of shocked Fe only in a small region further 
indicates the reverse shock reaches the Fe-rich ejecta very recently. These strong ejecta emission in the SE half is not 
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expected if the NW half has a higher surrounding ISM density and the SN exploded symmetrically in the center of the 
SNR. Therefore, it is very like ly that the SN explosion is also intrinsically asymmetric, as expec ted theoretically (e.g., 
Garcia-Senz fc WooslevI UMBh and also suggested in many previous observational works (e.g., iHamilton et al l [l997l: 
Uchida et al.ll2Q13[h 


5. SUMMARY AND CONCLUSIONS 

We develop new methods to conduct spatially resolved spectroscopy analysis in tessellated meshes each with total 
counts number above a given threshold. We apply this method to our XMM-Newton LP and archival XMM-Newton 
observations of SN1006, with total effective exposure times of 683, 710, and 439 ks of MOS-1, MOS-2, and PN. 
We extract spectra from 3596 tessellated regions each with a 0.3-8 keV total counts number > 10^. We fit these 
spectra with a simple “VNEI+SRCUT” model plus Gaussian lines representing O VII K6 — ( lines and sky background 
components. This simple spectral model with a 1-T thermal plasma is generally reliable to decompose the thermal 
and non-thermal emission and to characterize the average properties across the SNR. For the first time, we map out 
multiple physical parameters of the thermal plasma and the non-thermal component directly obtained fro m or derived 
based on the spectral fitting results. In particular, we adopt the shell-like geometric model introduced in iMiceli et al.l 
(j2Q12f ) and obtain an electron density map consistent with multi-wavelength estimates of the ambient ISM density. 

We also develop a new method to construct equivalent width (EW) maps based on the continuum interpolation 
with the spectral model of each tessellated regions. This method has the advantages of accounting for the curvature 
of the underlying continuum and decomposing the thermal and non-thermal components. We compare the EW maps 
constructed with this method and the EW maps constructed with a common linear interpolation method as well as the 
abundance map from our spatially resolved spectral analysis. The spatial distributions of emission lines as revealed 
by these three methods are qualitatively consistent with each other. 

We further extract spectra from some larger regions (than the tessellated meshes) and compute average parameters 
in them to confirm the prominent features as revealed by the parameter and EW maps. These features are generally 
preserved, although sometimes less prominent because of the mixture of the spectra extracted from different meshes. 
We analyze these spectra extracted from the larger regions with the same 1-T model as adopted for individual meshes. 
The spectral analysis results are generally consistent with the parameter maps, except for some “bad pixels” caused 
by the inadequate modeling of a few isolated meshes. 

In order to better characterize the average thermal and ionization states of SN1006, we construct probability distri¬ 
bution functions (PDFs) of kT and Uet of the thermal plasma. Previous measurements of kT and Uet in the literatures 
show large dispersions, but are all consistent with our measurements in different regions across the SNR. We thus con¬ 
clude that the diversity of the kT and Uet measurements in the literatures is at least partially caused by the different 
selections of spectral regions. We characterize the kT and Uet PDFs with several Gaussian distributions. Only with 
such a modeling could we characterize the average thermal and ionization states of such an extended source, the gas 
in which spans a large range of hydro dynamical evolutionary stages. 

For the first time, we present maps across the entire SNR of kT, Uet, Ton (ionization age with time dimension), Ue, 
and metal abundances of the thermal component, as well as a (radio-to-X-ray slope) and i'cutoff (cutoff frequency) 
of the non-thermal synchrotron emission, a is roughly constant in regions dominated by non-thermal emission, but 
i^cutoff shows clearly spatial variations and is higher in brighter non-thermal filaments. In general, the outer shells have 
low kT, Uet, and indicating that they are recently shocked ISM. The estimated electron number density represents 
a mixture of both shocked ISM and shocked SN ejecta, so could be adopted as upper limits of both components. The 
He map follows the broad-band X-ray surface brightness distribution, except for the two non-thermal limbs. There is 
a “Dark Belt” presenting on the Ue map, which also has higher Ne abundance, lower Uet, but comparable Uon than 
the surrounding regions. Based on the Ue map, we could estimate the total mass of the X-ray emitting plasma under 
the shell like geometric model, which is ^ 14 . 5/2 Mq. By comparing with the swept up ISM mass and the expected 
SN ejecta mass, we could further roughly constrain the volume filling factor / ^ 0.1. 

Spatial distribution of metals are traced with both the abundance maps and the EW maps constructed with either 
the linear EW method or the 2-D Spec EW method. O abundance is not significantly different from solar value over 
the entire remnant, and shows a low abundance hole roughly in the center of the SNR. The O VII K(5 — C/0 VIII EW 
ratio is significantly higher in this “O Hole” than any other regions. This “O Hole” and the surrounding regions are 
also the only part of the SNR with significant Fe emission, from both low-ionization Fe K lines and Fe L lines. All these 
features indicate that the “O Hole” should be an Fe-rich ejecta shell recently reversed shocked. However, the high Uet 
of this region also indicates that it should be largely comprised of shocked ISM. We therefore need a multi-temperature 
model to better decompose different thermal components. The “NW Shell”, which is known to be blast wave shocked 
high density ISM, shows significantly enhanced Ne emission (but the Ne abundance is still subsolar), but just moderate 
O and Mg emissions. We thus conclude that Ne is mostly ISM in origin, but O and Mg may be partly from the ejecta. 
Si and S emissions are stronger in the SE half of SN1006, and is clearly different from O, Ne, Mg, and Fe emissions 
in spatial distribution. The asymmetric distributions of metals strongly suggest that there is either an asymmetric 
explosion of the SN or the ISM distribution is highly asymmetric around SN1006. 
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APPENDIX 

ANALYSIS OF THE SKY BACKGROUND 

After removing all the bright point-like sources, we extract the spectra of sky background in an annulus with radii 
of 20' — 25' and centered at (ctj 2 ooo = 15^02"^54.03^, ^j 2 ooo = —41°55'39.86") (Fig. [T]). We fit the spectra with 
a “VMEKAL+power law” model subjected to the foreground absorption plus a second MEKAL model with solar 
abundance and no foreground absorption. Here the VMEKAL and MEKAL models represent thermal plasma in 
an equilibrium ionization state with or without relative abundances set. The first thermal component (VMEKAL) 
characterize the combined contribution from the possibly scattered light from the remnant as well as the background 
radiation from the Galactic corona. The abundances of some key elements, for example, O, Ne, Ee, are set free in order 
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to account for some prominent emission line features. The power law component, on the other hand, is used to char¬ 
acterize the possible contribution from the scattered li^h t of the non-therma l arcs and the cosmic background mainly 
consists of unresolved active galactic nuclei (AGN) fe.g.. iMoretti et ani2QQ3[) . The second thermal plasma component 
(MEKAL) represents the contribution from the local hot bubble, which is a structure in the solar neighbourhood, and 
is assum ed to have solar abundan ce and no foreground absorption. The temperature of this component is fixed at 
0.1 keV (jKuntz fc Snowdenll2QQQf) . Such a background analysis is not aimed at well describing the sky background 
physically, but at roughly characterizing different components of the sky background for subtractio n from the sou rce 
spectra. A similar double subtraction procedure is often used in the study of extended sources (e.g., iLi et al.ll200l) . 

Fitting the background spectra with the above model results in a temperature of the VMEKAL component of 
^ 0.2 keV an d subsolar abunda nces. These are consistent with the VMEKAL component of mainly Galactic halo 
in origin (e.g., iLumb et al.ll2002f ) and a negligible contribution from the scattered light of the SNR. The power law 
component has a photon index of ^ 1.4 —1.7 (slightly differen t for different instr uments), significantly smaller than that 
of the synchrotron emission from the non-thermal arcs (e.g., I Allen et al.ll2008[ ). but roughly consi stent with primarily 
cosmic in origin (from distant AGN) in both photon index and intensity (Fig. flSllLumb et al.ll2QQ2] ). This also indicates 
a negligible contribution from the scattered light of the non-thermal arcs of the SNR. 



Fig. 18.— Sky background spectrum (only PN to clarify) extracted from the whole annulus (black) and the four quadrants (colored) in 
Fig.[T] with point-like sources subtracted. The dotted lines show different model components of the whole annulus as described in the text, 
while the solid line is a combination of all these components. 


The above spectral analysis indicates a negligible contribution of the scattered light of the SNR to the sky background. 
The reason for this low s cattering is that the foreground H I (and so dust) column density toward SN1006 is low 
(Vh = 6.8 X 10^^ cm“^; iDubner et al.ll2QQ^ . given that SN1006 lies high above the Galactic plane. We further 
examine the azimuthal variation of the sky background to confirm this conclusion. X-ray emission of SN1006 shows 
significant azimuthal variation in both spectral shape and total intensity. Therefore, the sky background may have 
non-negligible azimuthal variations if the scattered light is important. We examine this possibility by extracting 
background spectra from each of the four quadrants of the annulus as shown in Fig. [U with azimuthal angles of 
13° — 103°, 103° — 193°, 193° — 283°, and 283° — 13° from west. As shown in Fig. [181 fhe flux variation in different 
quadrants is small, and could be largely attributed to the variation of the effective area scale of different regions. The 
scatter of the intrinsic intensity of the sky background is ^ 8% (at 1 a confidence level). Since the sky background only 
contributes a small fraction of the flux in the source regions (e.g.. Fig. n, this fluctuation in background flux should 
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have negligible effect on spectral modeling. Furthermore, we do not find any significant spectral variation in the four 
quadrants. In fact, the models of the four quadrants shown in Fig. [18] are exactly the same as that of the integrated 
spectra of the whole annulus, except for a constant normalization factor multiplied to all the model components. In 
conclusion, it is safe to apply the sky background from the whole annulus to all the source spectra. 


